Diverging Static Length Scales and Clustering in Cold Lennard-Jones Fluids 
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We report molecular-dynamics simulations on a three-dimensional two-component Lennard-Jones 
fluid. We identify two distinct static length scales. The longer length scale, which diverges at 
low temperatures, is continuous with the high-temperature static correlation length. The shorter 
length scale reflects static clusters that form below a T mc « 2 well above an apparent melting point 
T m ~ 1.08. We partially visualize the clusters; their maximum radius is temperature- independent. 
Cluster static and dynamic properties match Kivelson's model (S. A. Kivelson, et al., J. Chem. 
Phys. 101, 2391 (1994)) for frustration-limited cluster- forming glasses. 
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Many liquids, when cooled below their melting points, 
transform not into crystals but into amorphous glasses. 
There remains disagreement: Is glass formation an equi- 
librium thermodynamic transition or a dynamic effect? 
At a thermodynamic transition, static atomic position 
correlations would change markedly. Searches j4j,[f| of the 
radial distribution function g(r) of cold glass- forming liq- 
uids found no such changes. More recent studies 0,0 of 
multipoint dynamic correlation functions find increasing 
position correlations between particularly mobile or im- 
mobile atoms as the glass is approached. 

We performed molecular dynamics simulations on a 
three-dimensional Lennard-Jones fluid. In addition to 
the well-known local (r < 5 atomic diameters) order, we 
show that g(r) has extended (out to r > 11) order in cold 
systems. The crests and troughs of g(r) share a common 
envelope function; the function's length scale diverges as 
T — > 0. On extracting the extended order, cold Lennard- 
Jones fluids are found to form locally-ordered structures 
("clusters"), whose density is that of the bulk, and whose 
maximum radius is temperature-independent. 

Our fluid systems contained N = 125, 000 atoms, 
equally of species A and B, in an L x L x L (L « 43) box 
with periodic boundary conditions. The number density 
was p — NL~ 3 — 1.30. The potential energy was 
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with r the interparticle distance, i and j labelling species 



and Uij — for r > 2.5. Following Glotzer, et al. 1], we 
took o~aa — 1, o~bb — 5/6, and o~ab — 11/12. The inter- 
particle energy e and particle mass m were the same for 
all atoms. In natural units the temperature T is fce/e and 
the time unit is a ^(m/ 'e) 1 / 2 '. Numerical integration of 
Newton's equations of motion was accomplished via sym- 
plectic integration using the Calvo and Sanz-Serna fourth 
order methody, • The time step was never greater 
than 0.01, and less at T > 1. Simulations were made 
seriatim at 34 temperatures by progressive cooling and 
re-equilibration starting from an equilibrated T — 47800 
fluid. 



We calculated specific radial distribution functions 
9aa( t ), 9AB( r ), 9BB(f), and the total radial distribu- 
tion function g(r), with normalization g^- (oo) = 1, which 
give likelihoods of findingtwo atoms separated by r. Fol- 
lowing Steinhardt, et al.|8j, we also calculated spherical 
harmonic components 
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for atoms in a coordination shells, a coordination shell 
being a single peak of g(r) (region where g(r) > 1), with 
£li the vector from the center to atom i. To determine 
the const ant- volume melting point T m we heated crystals 
from T — 0.00. At T m , g(r) crosses over from solid to 
liquid behavior and atoms become mobile. The highest- 
melting-point crystal was body-centered-cubic, stable to 
T m « 1.08. 

Figure □ shows | g(r) - 1 | at T = 20, 2.01, and 0.56. 
Because we plot the magnitude, crests and troughs of 
g(r) all appear as peaks. At the high temperature T=20, 
order persists through four crests and matching troughs, 
with no order for r > 4. 

At T = 2, at least eight crests (and intervening 
troughs) appear in g(r). At T = 0.56 a dozen coordina- 
tion shells are visible; ordering extends to r « 11. This 
ordering is not microcrystallization: First, T m is 1.08, 
but ordering increases progressively starting at T > 10. 
Second, features of g(r) of warm crystals and the fluid do 
not match; for example, the second crest of g{r) is skewed 
toward small r in crystals, and toward large r in the 
fluid. Simulations on small (say, L < 20) systems might 
overlook extended ordering because periodic boundary 
conditions would overlap coordination shells on opposite 
sides of each atom and obscure their forms. Indeed, we 
dropped preliminary simulations with N — lb, 000 to en- 
sure we saw bulk properties. 

Most crests and troughs of g(r) fit accurately (but note 
Figure |5j) to 



g(r) =T 1 -T 2 (r- T 3 ) 2 exp(-T 4 r) 



(3) 
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FIG. 1: The all-atoms radial distribution function, plotted as | g(r) — 1 | against r, at temperatures 20(o), 2.01 (•), and 0.56(o). 
T = 20 and T = 2 data are displaced downward for clarity; | g(0) — 1 |= 1 always. Near-parabolas through T = 0.56 data 
are exemplary fits of eq.|3] omitted for the ninth crest and third and tenth troughs to permit judgement of data quality. Solid 
and dashed curves are envelopes based on eq. 0]and a simple exponential, respectively. The strong dependence of the range of 
liquid order on temperature is immediately transparent. 



Ti is the peak height, T% is the peak center, and T2 and 
T4 are shape parameters; figure^ shows exemplary fits. 

To characterize | g(f) — 1 | we fit peak maxima to en- 
velope functions, as done by Glotzer, et al. 1] with four- 
point two-time dynamic correlation functions. Our en- 
velopes were an exponential Tf t (r) = t Q exp(— nr) (where 
k is an inverse range) and 

rf t2 (r) = aexp(-6r c ), (4) 

with a, b, and c as fit parameters. The range of Tf* is 
defined to be k _1 or by volume integrals 

(R 2 ) = J 1 \ J . 5 

At large T, all peaks were fit. For T < 2, fits were 
restricted to peaks with r > 4. Crests and troughs of 
g(r) all fit to one envelope. 

The short-range (r < 4) part of g(r) of cold systems 
contains additional structure. Figure [21 shows deviation 
plots (peak height, minus envelope fitted for r > 4) at 
three temperatures. For r > 4 deviations are small. At 
T — 2.4 there are no deviations. Below T — 2, deviations 
increase markedly toward small r. At all T, deviations 
stop at the fourth coordination shell and do not continue 



to larger R. As a simple approximation we fit deviations 
to ~ exp(— K2r). The Hj 2 IS cL second length scale. (To 
reject artifacts: For T < 2 we also fit eq. 01 to all peaks, 
finding relative to the r > 4 fits a ten-fold larger error; 
the r < 4 deviations remained.) 

Figure [3] shows R, and k^ 1 against T. n~ l and 
R describe the outer coordination shells; they differ nu- 
merically because they weight large-r effects differently. 
k^ 1 describes inner coordination shells of cold systems. 
The increase in k^} corresponds to T- dependences of in- 
ner peak heights, all within the temperature- independent 
maximum cluster radius. These lengths scale as power 
laws 

t - T~ x . (6) 

for x — 0.30, 0.30, and 0.56, respectively. k _1 and R 
share an exponent; they reflect the same length scale. 
Exponents for k^ 1 and R differ; R and K2 represent phys- 
ically distinct lengths. We have thus found two distinct 
static lengths in cold Lennard- Jones fluids, both seem- 
ingly diverging at low temperature, contrary to any belief 
that static lengths are quiescent at small T. 

The scatter in Figure [2] at r < 4 is non-random. Fig- 
ure shows amplitudes T\ for the first peaks at vari- 
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FIG. 2: Deviations | g(r) — Tf t2 (r) | of peaks of g(r) at tem- 
peratures 0.56 (o), 1.0 (•), and 2.4, and (solid line) exp(— K2r) 
fit to r < 4.0 deviations at T = 0.56. 
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FIG. 3: Temperature dependence of (top to bottom) R, k~ 
and K% , and fits to £ ~ T~ x . 



ous T. Here T\ is normalized: First, division by T\ at 
T — 2 suppresses peak to peak height variations. Sec- 
ond, the T-dependence of the outer shells is suppressed 
by dividing out the T-dependence of the sixth crest of 
g(r), which is the most accurately known of the outer 
(r > 5) coordination shells. Below T — 2, the normal- 
ized Xi fall with decreasing T. Below T — 1.0, the first 
and second troughs and the second and third crests reach 
limiting amplitudes, but the first and second crests and 
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FIG. 4: Temperature dependence of the normalized peak am- 
plitude T\ (T, %) (i labelling the peaks of Figure |I] i = 1 being 
the excluded-volume core, and i = 12 being the sixth crest of 
g(r), for the early crests (labelled on abscissa) and interven- 
ing troughs (unlabelled on abscissa) of g(r), at temperatures 
(bottom to top) 2, 1.6, 1.2, 1.0, 0.86, 0.74, 0.69, 0.67, 0.65, 
and 0.55. Note inverted vertical scale. 



third and fourth troughs continue to decrease. It has 
long been known that the second peak has a character- 
istic glass-denoting characteristic, namely a notch and 
larger-r shoulder. We find weaker glass-denoting charac- 
teristics in other crests and troughs of g(r). 

Figure shows the detailed shapes of the inner crests 
and troughs of gAA(r) and ^(r) at our coldest tem- 
perature, 0.56. For both, the second crest (near r — 2) 
has a strong larger-r shoulder. Estimating the skewness 
of the peaks by the half- widths at half- height, the third 
trough is markedly skewed to smaller r, and the first 
and third crests are skewed to larger r. In contrast, at 
T — 2.01 (not shown), the first crest and trough are 
slightly skewed, but the outer crests and troughs are not 
skewed appreciably. 

We interpret the behavior of g(r) as reflecting two inde- 
pendent effects. First, there is a fundamental liquidlike 
order in g(r) that packs atoms in regular shells. As T 
declines, the regularity of the shells increases, leading to 
shells becoming visible at larger r. The range R of liq- 
uidlike ordering appears to increase indefinitely as T is 
reduced, reaching r!^llatr = 0.56. 

Second, below T = 2, short-range ordering appears. 
We interpret this effect as arising from the formation of 
local clusters. Clusters have a temperature-independent 
maximum extent, the fourth trough of g(r). Clusters 
are partially visualized in Figure El (T dependences of T\ 
of crests one and two and trough three), and Figure 
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FIG. 5: Expanded view of g(r) (lines), qaa (open points), 
and cjbb filled points) at T = 0.56, showing emergent line 
structure. g(r) omitted for the second peak and third trough. 
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FIG. 6: Dynamic correlation (Y2K) between a spherical har- 
monic component of the local density of the first coordinate 
shell and the kinetic energy at temperatures 0.59(+), 0.80 
(solid line), 1.20(dashed line) and 2.00(o). 

(shoulders and asymmetries in the first three crests and 
third trough of g(r) and the ga(r)). 

Third, clusters are transparent to long- wavelength 
scattering. At T = 2 and T = 0.56, J Q 4 (^(r) - l)dr » 0. 
Cluster formation rearranges atomic positions for r < 4 



but do not change the cluster mass. 

Within a sphere of diameter 8, a cluster might have 300 
atoms. However, all atoms in a sphere need not be part of 
a cluster. Donati, et al.|6j concluded from a four-point, 
two-time distribution function that the radial distribu- 
tion function for pairs of immobile atoms, in a system 
and temperature close to ours, is relatively well-ordered 
at small distances and fades to zero beyond r « 4, the 
size of our clusters. Our static clusters are plausibly the 
immobile dynamic clusters of Donati, et al. 

Observed clusters have properties predicted by Kivel- 
son, et al.[{| for frustration- limited cluster- forming 
glasses, namely: Clusters (i) are stable at temperatures 
above the melting temperature, (ii) manifest frustration 
by having a temperature-independent maximum size, 
and (iii) have a local order orthogonal to the density. Ap- 
plication of the Kivelson glass model is complicated by 
the low-T divergence of the liquidlike length scale, and 
by T-dependent cluster internal structure seen in K2- 

Finally, Kivelson, et al. predict that glass forma- 
tion creates dynamic correlations between local order 
parameters orthogonal to the density, and time deriva- 
tives A of principle mechanical variables. We evaluated 
(Y2K) — (Q2o(0)K(t)), where Q2o(0) was computed for 
the first coordination shell, and K(t) — p 2 /2m is the ki- 
netic energy of the central atom at later time r. Figure 
©shows that at small T (Y2K) gains (inset) new r ^ 
structure, as predicted by Kivelson, et al.'s model 9]. 
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